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<^ , Abstract 

We present a new Monte-Carlo method for estimating the chemical potential of model 
polymer systems. The method is based upon the gradual insertion of a penetrable 'ghost' 
polymer into the system and is effective for large chain lengths and at high densities. In- 
sertion of the ghost chain is facilitated by use of an expanded ensemble in which weighted 
' transitions are permitted between states characterising the strength of the excluded vol- 

ume and thermal interactions experienced by the ghost chain. We discuss the imple- 
mentation and optimisation of the method within the framework of the bond fluctuation 
model, and demonstrate its precision by a calculation of the finite-size corrections to the 
OS ■ chemical potential. 
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1 Introduction 



Measurements of chemical potential by computer simulation are an important tool for deter- 
mining the phase behaviour of model systems. For many years the standard technique for 
achieving this in atomic and simple molecular systems, has been the Widom particle insertion 
method []J . This method involves repeatedly sampling the potential Ut experienced by a single 
particle when it is added at random to the system of interest. The excess chemical potential 
\i ex is then given by 



PfMea = - In < exp(-pU t ) > (1) 

For polymeric systems, however, the basic Widom method cannot be employed for obtaining 
estimates of fi ex . The extended structure of polymers renders the method unsuitable for all 
but the shortest chains and lowest densities, there being a prohibitively small probability 
of inserting a polymer at random without violating excluded volume restrictions. Clearly 
therefore different approaches are called for, and indeed several new techniques for measuring 
the chemical potential of model polymer systems have recently been developed ||. We describe 
them in turn. 

A method that ameliorates some of the difficulties associated with random chain insertion 
is the so-called configurational bias Monte Carlo method, first proposed by Siepmann || and 
later generalised by Frenkel and coworkers and de Pablo and coworkers |J. This technique 
utilises a biased insertion technique to 'grow' a polymer of favourable energy into a frozen 
snapshot of the system, thereby extending the range over which the Widom insertion technique 
is viable. A book-keeping scheme maintains a record of the statistical bias associated with 
choosing energetically favourable chain conformations, and this bias is subsequently removed 
when the chemical potential is calculated. The method is also amenable for use with Gibbs 
ensemble simulations, facilitating direct study of polymer phase equilibria |], [7| . However in 
practice, the utility of the configurational bias approach is limited since it tends to break down 
for longer chain lengths (> 20 monomers) and higher densities (p > 0.6). 

Another technique for obtaining fi ex is the Chain Increment method of Kumar and cowork- 
ers 0, H ||. This method hinges on the observation that the chemical potential of a chain 
of length N is equivalent to the sequential sum of the incremental chemical potentials of all 
successively shorter chains N — 1, N — 2, .... Since the incremental chemical potential seems 
to become largely independent of chain length for longer chains (N > 15) JT0], it is possi- 
ble to estimate the total chemical potential for a chain of arbitrary length by calculating the 
incremental chemical potential associated with adding one monomer to a chain end, and multi- 
plying by the chain length. Separate sets of simulations of short chains (N < 15) may then be 
performed to provide corrections associated with the small- N dependence of the incremental 
chemical potential. This method thus allows estimates of the chemical potential for arbitrary 
chain lengths, and indeed has facilitated the calculation of the liquid-vapour phase diagram 
for a bead-spring polymer model |TT] . 



Recently Miiller and Paul have developed a novel approach to calculating the chemical 



potential of polymer chains [12]. Their method involves a thermodynamic integration over 
the excluded volume interaction of a penetrable 'ghost' chain immersed in a system of 'real' 
polymers. A control parameter A tunes the strength of this interaction such that for A = 
1 the ghost decouples from the normal polymers, while for A = it behaves as a normal 
polymer. Separate simulations are performed at each of several A values in the range (0,1), 
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and thermodynamic integration with respect to A, aided by multi-histogram extrapolation 
techniques, are used to obtain the chemical potential. The method was successful in obtaining 
p ex for lattice polymers up to iV = 80 monomers and at melt densities (p = 0.5). 

In this paper we present a further method for calculating the chemical potential of poly- 
mer chains that is effective both at large chain lengths and high densities. The method is 
efficient, requiring only a single simulation, and exact in the sense that no approximations are 
involved. In the following sections we describe the method and discuss its implementation and 
optimisation. 



2 Method 
2.1 Background 



The strategy of the present method is similar in spirit to that of Muller and Paul JT2|. We 
consider a system of volume V = L d containing n polymers of length N, plus a single penetrable 
ghost polymer. The monomers of the normal polymers experience an attractive short-range 
interaction J(r) and are subject to strict excluded volume restrictions. An intra-molecular 
bond potential U(r) controls the bond lengths and angles between neighbouring monomers of 
the chain. 

The ghost polymer is subject to the same intra-molecular bond potential as the normal 
polymers and its monomers must not overlap with one another. However, the excluded volume 
restriction between ghost monomers and normal monomers is now replaced by an athermal 
effective repulsive potential K = —\ln\, where A is a parameter, variable in the range (0, 1), 
that controls the strength of the repulsion. The magnitude of the attractive thermal interaction 
between ghost monomers and normal monomers is also modified by a factor (1 — A). 

For a given choice of A, the canonical partition function of the system is given by 

(n+l)N 

Z{n,N,V,\)= J] { / drA e-^ Et ^ +E ^ (2) 

where the integration is over all configurations of the monomers i, subject to excluded volume 
restrictions. 

The thermal part of the configurational energy E t ({r}) is defined by 



Et({r}) = E Jdti ~ L k \)8 aj0lk + (1 - A) £ J(\ Lj - r fc |)(l - 5 ajCtk ) + £ Ufa-rJ (3) 

(j,k) (i,fe) (I,m) 

where j and k run over all monomers, while I and m run over all neighbouring pairs of monomers 
on the same chain. We assign aj = 1 for a normal monomer and ctj = for a ghost monomer. 
The athermal part of the configurational energy E a ({r}) is given by 

E a ({r}) = -- p N {{r_})ln\ (4) 

where N ({r}) is the total number of overlaps between the 'excluded volume' shell of ghost 
monomers and those of normal monomers. 

Clearly, by tuning A one is able to interpolate smoothly between two special cases: 
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• A = 0, for which the ghost behaves as a normal polymer 

• A = 1, for which the ghost decouples from the polymer system. 

The ratio of the partition function at these limits defines the excess chemical potential: 



Following the recently proposed method of expanded ensembles [fT3| , p3| , we now consider an 
ensemble that facilitates direct measurement of the ratio of the partition functions appearing 
in equation [|[ The ensemble is defined by the partition function: 

M-l 

tt = Z m (n, N, V, A m ) exp(?7 m ) (6) 

m=0 

where < m < (M— 1) indexes a set of M A-states spanning the range (0, 1), with Ao = and 
Xm-i — 1- The Z m are the canonical partition functions defined in equation ^ with A = A m , 
while the r] m are positive pre-weighting factors the significance of which will be described. 

To realise a Monte- Carlo simulation within this expanded ensemble, it is necessary to allow 
for two types of MC steps. Firstly there are the usual trial monomer displacements at constant 
A. Secondly there are trial transitions between neighbouring A-states, the spatial coordinates 
of all polymers remaining fixed. These A-transitions are accepted or rejected according to some 
criterion satisfying detailed balance, such as metropolis: 

P a (n -> m) = min(l,exp[N (hiX m - lnA n ) + (3{X n - \ m )E G +r] m - r] n )) (7) 

where E G = J2(j,k) J{\Lj — £fc|)(l — ^cj is the interaction energy of the ghost polymer with 
the other polymers. It should be stressed, however, that these A-transitions do not drive the 
system out of equilibrium as they would if we were considering a simple canonical ensemble. 
Within the expanded ensemble framework, both types of moves — monomer displacements and 
A-transitions — are necessary to bring the system into equilibrium. 

In the course of the expanded ensemble simulation a histogram p(X m ) can be accumu- 
lated describing the probability with which each of the M A-states is visited. Formally this 
probability distribution takes the form: 



From which it follows that 



p(X Q ) Z exp(?7o 

The excess chemical potential is then calculable as: 
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p(Xm) = -^exp(?7 m ) (8) 



p(X M -i) _ Z M _x exp(77M-i) ^ 



Hex = In 

where 



P(X 



M-l 



(10) 



P(Xo) 

P(X m ) = p{X m ) exp(-r? m ) (11 
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Now, the principal difficulty in implementing such a scheme, is that the ratio of the partition 
functions Zm-i/Zq appearing in equation |^ can (even for short chains) span very many orders 
of magnitude. In order to obtain a histogram p(X m ) having counts for all m, it is therefore 
necessary to employ a pre-weighting scheme that will encourage the simulation to sample A- 
states that would otherwise be highly unfavourable. This is achievable by suitable choices of 
the pre-weighting factors 7] m featuring in equations |6|-[TT|. Unfortunately however, the choice 
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of suitable rj m cannot be determined directly a-priori. For example, the weights that would 
yield a perfectly flat sampled distribution (i.e. p(X m ) = constant), are given by 

f 1 — In Z(n,N,V,X)dX (12) 

in which Z(n, N, V, A) is, of course, unknown. Nevertheless, it turns out that it is possible to 
implement a straightforward iterative scheme for obtaining suitable r\ m . 

2.2 Determining the pre-weighting factors 

The goal is to determine values for the pre-weighting factors r] m that allow the system to 
visit each A-state with probabilities of similar magnitude. This may be accomplished by the 
following procedure. 

1. Choose a set of M A values at which states are to be placed. The exact value of M is 
not too important at this stage, although it should not be too small. A discussion of the 



optimal choice for M is given in section 3.1. For most purposes, however, between 5 and 



15 states will suffice. For reasons that will also be described later, these states should be 
distributed such that their density increases with increasing A. 

2. Starting with all pre-weighting factors 7] m = 0, conduct a relatively short simulation 
(10 4 -10 5 Monte-Carlo sweeps), accumulating the histogram p(X m ). As a rule, the vast 
majority of counts will be confined to the end state m — M — 1. Provided, however, that 
the adjacent state (m = M — 2) is placed close to A = 1.0, the system will also (albeit 
rarely) make some excursions to this state too. 

3. Suppose now that counts were obtained for p(X m ) in the i states M — i, M—i+1, .., M — 1. 
New estimates for r] m are obtained by linear extrapolation using the iterative formula: 



V' m = r] m -\np(X m ) M -i<m<M -1 (13) 

V'm= Vm—i + [fe^gj Hu-i - v' M -i + i) 0<m<M-t 

4. Next perform another short simulation with the new pre-weighting factors r)' m . In this 
second simulation, p(X m ) will receive counts in at least i + 1 A-states. 

Stages 3 and 4 are then simply repeated until weights are obtained that allow the simulation 
to explore the entire A domain. This iteration procedure can of course be easily automated 
if desired. Once suitable weight factors are obtained, a longer simulation can be performed 
to obtain the desired statistical accuracy for p(X m ). The chemical potential is then obtained 



directly from equation 10 . 
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In some circumstances shortcuts to obtaining suitable weight factors can also be exploited. 
For example, it is advantageous to carry out the above procedure for a small system. This 
entails little computational effort, and assuming that the finite-size dependence of [i ex is not 
too great, the weights acquired should serve for long runs on larger systems. Additionally, if it 
is desired to obtain \i ex for a variety of chain lengths at a single density, weights determined for 
a rather short chain can be extrapolated to longer chain lengths by exploiting the approximate 



linear dependence of fi ex (and r] m ) on the chain length fllP] . Refinements to the extrapolation 
procedure outlined above can also be envisaged. For instance, a polynomial fit to all the 
weights already determined, might extend the range of the extrapolation procedure, allowing 
two or more weights to be determined per iterative cycle instead of just one. In practice, 
however, (and as is demonstrated in the following section), the above procedure yields the 
weight factors surprisingly rapidly. 



3 Application to the Bond Fluctuation Model 

In order to test our method for measuring fi ex , we have implemented an expanded ensemble 
simulation using the bond fluctuation model. This is a lattice model for polymers in which 
polymer chains are represented as N monomer units, each occupying a unit-cell of a simple 
cubic lattice. Successive monomers are connected by bonds, the bond vector of which is allowed 
to assume one of 108 distinct values. The number of allowed bonds is dictated by prescribed 
restrictions on the bond lengths and the necessity to prevent chains from crossing one another. 
Aside from these restrictions however, no explicit intra-molecular potential is incorporated. 
Excluded volume effects are catered for by requiring mutual self avoidance of monomers: no 
lattice site can be occupied simultaneously by two monomers (unless, as is possible in the 
present case, one happens to be a ghost monomer). For simulations in the thermal regime, an 
attractive inter-molecular potential can be applied as required. Polymer moves are facilitated 
either by a local monomer displacement algorithm (as was the case in the present work), or 
alternatively by reptation type moves. Further details concerning the model can be found in 
references KT5I, [IB . 



Simulations were performed using a metropolis algorithm for both the local monomer dis- 
placement algorithm and the A-transitions |T7|. In the course of the simulations, a number of 
different chain lengths and densities were studied in both the thermal and athermal regimes. 
We begin, however, by describing how the weight factors r] m were obtained for an athermal 
system of polymers of length N = AO monomers in a volume V = 40 3 at volume fraction 
p = 0.4. 

A total of eight A-states were employed, these being placed at the A values shown in table [H 
It is not possible to place a A-state exactly at A = since this causes a computational error 
when the transition probabilities (equation [1]) are calculated. Nevertheless, the value of Aq can 



be made extremely small, resulting in negligible error. For reasons made clear in subsection [3[T 
the distribution of A-states was chosen such that their density increased with increasing A . 
The ratio of the number of displacement moves per monomer to the number of attempted 
A-transitions was set at 4 : 1 respectively. 

The evolution of the iterative procedure by which the weights were obtained, is also shown 
in table |l|. Each run constitutes one iterative cycle, in which a further weight factor was 



determined. It is seen that for each run the extrapolated guess (equation [H5p systematically 
overestimates the r) m value to be determined, resulting in proportionally more weight for the 
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corresponding entry in p(X m ) on the ensuing run. However, this overestimate is subsequently 
corrected on the following cycle, when the weight factor is reduced with respect to the others. 
The final weights are presented in the last column of table [TJ. Despite the need for several 
cycles to determine the weight factors, the whole procedure was surprisingly rapid, consuming 
only 30 minutes CPU time on an IBM RS/6000 workstation. 

Having obtained suitable weight factors, a longer run of 2 x 10 6 Monte Carlo sweeps (MCS) 
was performed. Figure |l](a) shows the pre-weighted distribution p(X m ) obtained from this run, 
while figure [L](b) shows the form of the reweighted distribution p(A m ) = exp(— r] m )p(X m ). The 
associated chemical potential is p ex = 61.4(1), corresponding to a random insertion probability 
of 10- 26 . 

With regard to figure 0(a), it should be pointed out that the apparent noise in p(X m ) is 
illusory. The variations in p(X m ) represent not statistical uncertainty (which in this instance 
is much smaller than the variation), but rather deviations of the weights from their 'ideal' 
values as prescribed by equation [T2|. Provided though that these deviations are not so great 
as to hinder the random walk of the system between the A-states, there will be practically no 
influence of the weight factors themselves on the statistical accuracy of p ex , since their effect 
is simply removed when the reweighted distribution p(X m ) is formed. 

To demonstrate that the method is effective for even longer chain lengths and higher den- 
sities, we have studied an athermal system of polymers of chain length N = 80 at a volume 
fraction p = 0.5, contained in a volume L = 40 3 . The weight factors were determined exactly 



as before, although as discussed in subsection |3Tl~| , more A-states (14 in total) were required 



to prevent the acceptance rate falling too low. The ratio of displacement moves per monomer 
to attempted A-transitions was set at 10 : 1 respectively. Having determined suitable weight 
factors, a production run of 4 x 10 6 MCS was carried out, resulting in the distributions shown 
in figure |2|. The corresponding estimate for the chemical potential is [i ex = 175.9(1), in good 
agreement with that obtained previously under the same conditions using thermodynamic 



integration fl2 | 



Turning now to the thermal regime, we have employed a short-range square well potential 
to study a system of volume L = 40 3 and chain length N = 40 at a volume fraction p = 0.2 
and temperature T = 2.0. The procedure for determining the rj m is similar to before, although 
a complication arises in the fact that p(X m ) is no longer monotonic. The forms of p(X m ) and 
p(X m ) are shown in figure [5[ The latter distribution displays a broad minimum, manifesting 
the competition between the positive athermal contribution to the chemical potential (which 
dominates as A — > 1) and the negative thermal contribution (which dominates for small A). To 
obtain the weight factors in this case it was necessary to start the extrapolation procedure at 
both ends of the A range, (i.e. m = and m = M — 1) and work in towards the position of the 
minimum. The final complete set of weights was then obtained by requiring continuity at the 
A value where the two sets coincide (in effect by adding a constant to one set). The chemical 
potential estimated from a run of 2 x 10 6 MCS is p ex = —0.01(7). This very small value of p ex 
indicates the close proximity of the temperature to the O-point of the model. 



3.1 Optimisation of the method 

In situations where very high statistical accuracy is required, long runs can be necessary and a 
little effort to optimise the method may pay dividends. Two interconnected factors influence 
the statistical accuracy that can be realised for a given expenditure of computational effort. 
Firstly there is the choice of the number of A-states M, and secondly there is the acceptance 
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rate for A-transitions. These factors jointly control the magnitude of the correlation time for 
the A-transitions. To the extent that p(A TO ) can be regarded as flat, the system will execute a 
one-dimensional random walk between the A-states and the correlation time will be given by: 



r(M) 



M 2 



(14) 



Pa(M) 



where r(M) is measured in units of transition attempts and P a (M) is the acceptance rate for A 
transitions, which we assume to be constant for each state. To appreciate how P a (M) depends 
on the number of A-states M (and their placing) it is instructive to consider the transition 
probabilities for an athermal system: 



In addition to ensuring that p(X m ) is approximately flat, the role of the weights in this equation 
is to partially compensate for positive changes in the exponent associated with transitions to 
smaller A values. The extent to which they are successful in this regard (and hence the size 
of the acceptance rate) depends on the degree of overlap between the distributions of iV a (A m ) 
and iV (A n ). In order to ensure a sufficiently high acceptance rate, sufficient A-states must 
be employed to guarantee an appreciable overlap between the distributions. Moreover, to 
maintain an approximately constant acceptance rate for all states it is necessary to increase 
their density at higher A values. This is illustrated in figure |3] where the distribution function 
N (X) is plotted for a variety of fixed A values. The distributions shown have approximately 
equal overlap although the A-states are placed closer together for larger A. 

The above considerations illustrate the compromise to be struck in minimising the corre- 
lation time: if M is too large, r will be large and accordingly it will be difficult to accumulate 
uncorrelated statistics. Similarly, if M is too small, a low acceptance rate will ensue and con- 
sequently r will again be large. However, it should be stressed that the correlation time is not 
the sole factor governing the optimum choice of M . For a given length of run, the statistical 
quality of p(A m ) may be improved by choosing a smaller value of M, even at the expense of 
a larger correlation time. If the sampling frequency is chosen commensurate with the correla- 
tion time (equation |13|), a fraction 2/M of the total number of samples will fall into the two 
extrema states p(A ) andp(AA/_i) used to calculate \i ex . Now, since the statistical error on \i ex 
is dependent on the number of entries in these end states, it follows that the optimal choice of 
M is given, not by the minimum in r(M), but by the minimum in M 3 /P a (M). In figure [5] we 
plot this function for a system of athermal polymers with iV = 40 and p = 0.2. A minimum 
obtains for M — 6, corresponding to an acceptance rate P a ~ 20%. Thus it would appear that 
in the present method, the optimal acceptance rate can be considerably less than the value of 
50% or so used for normal particle displacements. 

Finally we address briefly the question of the optimal choice of the ratio of local displace- 
ment moves to A-transitions. Clearly the issue is again one of efficiency. If too few local moves 
are attempted between transition attempts, then the number of overlaps N Q will systematically 
deviate from the canonical average < N(X) > and the acceptance rate may suffer. On the 
other hand, if the ratio of local moves to transition attempts is too high, a longer run will be 
necessary to obtain good statistics in p(A m ). For the present model, we find that following a 
A transition, between 4 (N = 40) and 10 (N = 80) local moves attempts per monomer are 
sufficient to decorrelate the number of overlaps. We note also in this context, that relatively 
few local moves are required to decorrelate the spatial position of the ghost polymer. This is 



P a (n -»• m) = min ( 1, exp [N ( In A m - lnA n ) + r] m - r) n ]) 



(15) 
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because the ghost polymer frequently samples A-states in which it completely decouples from 
the normal polymers, thereby allowing it to diffuse around the system with relative ease. This 
latter feature of the method is in favourable contrast to the thermodynamic integration method 



12], in which several hundred trial displacements are required for a ghost polymer with fixed 



small A to diffuse away from its starting point. 

3.2 Finite-size corrections to the chemical potential 

Measurements of the chemical potential in finite-size systems are known to be relatively 
strongly effected by finite-size effects. Indeed for a system of n particles, the leading cor- 
rections vary like n _1 . Recently, Siepmann et al have derived an expression for the leading 
finite-size dependence of p ex in terms of density derivatives of the pressure |TH 



A/i ex .(n) 





k h T 777^ - phT 



(d 2 P/dp 2 ) 
(dP/dp) 2 



0(n' 2 ) (16) 



where p = n/V is the polymer number density. Provided that the equation of state is known 
independently, this equation permits a calculation of the finite-size dependence of the chemical 
potential, thereby facilitating an extrapolation to the thermodynamic limit. 

As a test of the precision and accuracy of our method, we have attempted to measure the 
finite-size dependence of p ex for our model and compare it with equation [16|. To this end we 
have simulated systems of athermal polymer chains of length N = 40 at volume fraction p = 0.2 
for system sizes L = 40, 47 and 60. In accord with the findings of the previous subsection, 6 
A-states were utilised. Runs comprising 10 7 MCS were performed, and the A value was sampled 
every 180 transition attempts. For each system size studied, 16 independent runs were carried 
out in order to test the statistical independence of the data and assign errors to the results. 

The resulting estimates for the excess chemical potential are // ex (40) = 19.78(2), /x ea .(47) = 
19.69(2), p ex (60) = 19.59(2), where the value for L = 47 has been corrected for the fact that 
p = 0.2 does not correspond to a whole number of polymers. Figure shows how these 



estimates compare with equation [16]. The broken line represents the predicted finite-size shift 
in p ex obtained by feeding the equation of state data of Miiller and Paul [0 into equation |T6"| . 
Assuming an infinite-volume chemical potential of p ex (oo) = 19.50, the measured finite-size 



shift in p ex is consistent to within error with the prediction of equation |T6. 



4 Conclusions 

In summary, we have presented an efficient method for measuring the chemical potential of 
polymeric systems that is effective at both large chain lengths and high densities. Although 
some preliminary effort is required to find suitable pre-weighting factors, this is not excessively 
time consuming and is generously rewarded by the very high accuracy that can subsequently 
be attained. In contrast to the Chain Increment method M, the estimate of p ex derives from 
a single simulation and is in principle exact. The present method should therefore provide a 
more accurate alternative in applications such as determining the polymer- solvent coexistence 
curve of long chain molecules (where configurational bias Gibbs ensemble simulations are not 
feasible). Finally we remark that the present approach is not limited to homopolymer chains, 
and should also facilitate measurements of p ex for copolymers as well as more complex polymer 
architectures such as branched or ring structures. 
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Figure 1: (a) The pre- weighted distribution p(X m ) for a system of athermal polymers of chain 
length N = 40 and volume fraction p = 0.4, contained in a volume V = 40 3 . (b) The reweighted 
distribution p(X m ) = exp(—r] m )p(\ m ) expressed on a logarithmic scale. 



Figure 2: (a) The pre-weighted distribution p(\ m ) for a system of athermal polymers of chain 
length N = 80 and volume fraction p = 0.5, contained in a volume V = 40 3 . (b) The reweighted 
distribution p(X m ) = Gxp(—r] m )p(\ m ) expressed on a logarithmic scale. 



Figure 3: (a) The pre-weighted distribution p(X m ) for a system of thermal polymers of chain 
length iV = 40, p = 0.2 at temperature T = 2.0, contained in a volume V = 40 3 . (b) The 
reweighted distribution p(X m ) = exp(—r] m )p(\ m ) expressed on a logarithmic scale. 
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Figure 4: Distributions of the number of overlaps between ghost monomers and normal 
monomers at various fixed A. 



Figure 5: The measured form of the function M 3 /P a (M) 



Figure 6: The measured value of the chemical potential for a system of athermal polymers 
of chain length N = 40 at volume fraction p = 0.2 for system sizes L = 40,47,60. A con- 
stant corresponding to the infinite volume estimate p ex = 19.50 has been subtracted from the 
data. Also shown (broken curve) is the predicted finite-size dependence of p ex following from 
equation [16] utilising predetermined equation of state data [O . 





Run 1 


Run 2 


Run 3 


Am 




p(A m ) 




p(A m ) 




p(A m ) 


A = 0.00001 


0.00 


0.00 


12.60 


0.000 


24.28 


0.000 


Ai = 0.35 


0.00 


0.00 


12.60 


0.000 


24.28 


0.000 


A 2 = 0.55 


0.00 


0.00 


12.60 


0.000 


24.28 


0.000 


A 3 = 0.7 


0.00 


0.00 


12.60 


0.000 


24.28 


0.000 


A 4 = 0.8 


0.00 


0.00 


12.60 


0.000 


24.28 


0.567 


A 5 = 0.9 


0.00 


0.00 


12.60 


0.423 


12.18 


0.168 


A 6 = 0.97 


0.00 


0.022 


3.79 


0.299 


3.72 


0.142 


A 7 = 1.0 


0.00 


0.977 


0.00 


0.277 


0.00 


0.113 





Run 4 


Run 5 


Run 6 




Am 




p(A m ) 


Vm 


p(A m ) 




p(A m ) 




A = 0.00001 


33.63 


0.00 


47.13 


0.00 


57.78 


0.009 


61.02 


Ai = 0.35 


33.63 


0.00 


47.13 


0.00 


57.78 


0.673 


53.67 


A 2 = 0.55 


33.63 


0.00 


47.13 


0.656 


45.03 


0.060 


43.90 


As = 0.7 


33.63 


0.423 


32.77 


0.052 


33.21 


0.063 


32.50 


A 4 = 0.8 


22.75 


0.113 


23.21 


0.064 


23.44 


0.042 


23.14 


A 5 = 0.9 


11.86 


0.133 


12.16 


0.066 


12.37 


0.036 


12.21 


A 6 = 0.97 


3.57 


0.153 


3.73 


0.079 


3.75 


0.034 


3.64 


A 7 = 1.0 


0.00 


0.177 


0.00 


0.082 


0.00 


0.031 


0.00 



Table 1: The weight determination procedure described in the text for N=40, p 
runs are 50000 MCS. 



0.2. All 
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